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ABSTRACT 


This  thesis  develops  the  optimal  control  law  necessary  for 
a minimum  time  escape  from  a high  Earth  orbit  utilizing  solar 
radiation  pressure.  A force  model  Is  developed  for  an  Ideal 
flat  sail  and  a proposed  Jet  Propulsion  Laboratory  sq\iare  sail. 

A two-body  Inverse  square  force  field  model  Is  assumed.  The 
Earth's  orbit  about  the  Sun  Is  assumed  clrculeu*»  and  the  solar 
flux  Is  assumed  constant  during  the  escape  maneuver.  The  Initial 
state  Is  given,  and  the  only  condition  on  the  final  state  Is 
that  the  energy  of  the  spacecraft  equals  zero.  A modified 
Newton-Raphson  method  Is  used  to  solve  the  two-point  boundary 
value  problem.  It  Is  found  that  the  performance  of  any  solar 
sail  escape  trajectory  tends  to  fall  between  the  limiting  cases 
of  similar  polar  and  ecliptic  escape  trajectories.  A velocity- 
dependent  control  law  Is  also  examined  and  found  to  be  a good 
approximation  to  the  optimal  control  law. 
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1.  INTRODUCTION 


1 

( 


1.1  BackPTOiind  j 

A solar  sail  is  a device  which  creates  a propulsive  force 

from  solar  radiation  pressure.  The  solar  radiation  pressure 
results  from  changes  In  the  momentum  of  Incident  photons  on  the 
sail's  reflective  surface.  This  force  can  be  used  to  Increase 
or  decrease  the  energy  of  the  solar  sail*  thereby  allowing  travel 
In  the  vicinity  of  the  radiation  source. 

One  of  the  first  examinations  of  a solar  sail  powered 
vehicle  supposedly  occurred  In  a paper  by  the  Russian  author, 

F.A.  Tsander,  In  1924(ref  3).  It  was  not  until  the  late  1950* s 
that  serious  Interest  In  the  concept  of  solar  sail  propulsion 
occurred.  The  majority  of  the  papers  written  on  solar  sail 
propulsion  have  been  concerned  with  the  solar  sail  as  a means 
of  Interplanetary  travel.  The  use  of  a solar  sail  as  a rneauis 
of  planetary  escape  has  been  considered (ref  4 and  11),  but  the 
analysis  has  been  rather  restricted. 

Since  the  Space  Shuttle  program  Is  envisioned  to  make  low 
Earth  orbits  readily  attainable,  the  solar  sail  Is  at  present 
being  considered  for  various  high  Earth  orbit  and  Interplanetary 
missions.  The  solar  sail  Is  rather  Interesting  In  regards  to 
Interplanetary  travel  In  that  It  has  the  potential  for  Improved 
performance  over  Interplanetary  ballistic  and  solar  electric 
propulsion  systems.  Inherent  to  any  Interplanetary  mission  Is 
the  successful  escape  from  Earth. 

1.2  Ob.lectlve 


source,  the  optimal  escape  trajectory  from  Earth  should  be  the 
trajectory  that  reaches  escape  In  the  minimum  time,  A straight 
numerical  Integration  of  the  equations  of  motion,  in  order  to 
determine  the  optimal  trajectory,  usually  results  In  numerical 
problems.  These  numerical  problems  are  due  to  a large  escape 
time (In  the  order  of  months)  and  to  the  small  Integration  time 
step  needed.  The  method  of  averaging  the  orbital  elements  as 
discussed  In  references  8, 9 tend  10,  has  proved  successful  In  the 
determination  of  optimal  orbit  raising  trajectories.  As  the  ratio 
of  the  propulsive  force  to  the  local  gravitational  force 
becomes  large,  the  method  of  averaging  becomes  Inaccurate,  For 
the  solar  sails  discussed  In  Appendix  A,  It  appears  that  the 
method  of  averaging  should  be  terminated  when  the  trajectory 
reaches  the  vicinity  of  100000  kilometers  from  the  center  of  the 
Earth,  It  Is  therefore  the  purpose  of  this  thesis  to  develop 
an  analytical  model  that  will  calculate  the  control  law  necessary 
for  a minimum  time  escape  from  a high  Earth  orbit  with  a solar 
sail.  It  Is  felt  that  a model  similar  to  the  one  presented  In 
this  thesis  could  be  combined  with  an  averaging  method  to  deter- 
mine the  complete  optimal  escape  trajectory  of  a solsur  sail 
from  a low  Earth  orbit. 

Some  basic  ass^lmptlons  are  made  In  the  solution  of  this 
problem* 

1.  The  state  Is  completely  given  at  the  time  of  launch. 

2.  The  solar  flux  Is  constant  during  the  escape  maneuver, 

3.  The  distance  between  the  Sun  and  the  spacecraft  Is  equal 
to  the  distance  between  the  Sun  and  Earth. 
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4,  The  Earth's  orbit  about  the  Sun  Is  circular, 

5,  The  Earth  has  an  Inverse  square  gravitational  field, 

6,  The  only  gravitational  field  which  affects  the  space- 
craft is  that  of  the  Earth, 

7,  The  effects  of  Earth  shadowing  is  neglected, 

1,3  Synopsis 

This  thesis  is  divided  into  five  chapters.  Chapter  one  is 
the  introduction.  Chapter  two  developes  an  equation  for  the 
magnitude  and  direction  of  the  resultant  force  on  a solar  sail 
due  to  solar  radiation  pressure.  Chapter  three  developes  the 
differential  equations  of  motion  for  the  state  and  co-state. 
Optimal  control  theory  is  applied  to  determine  expressions  for 
the  time  history  of  the  control  variables  and  the  necessary 
boundary  conditions  on  the  co-state.  Chapter  four  explains  some 
of  the  results  of  optimal  control  escape  trajectories.  Two 
other  types  of  control  laws,  open-loop  control  and  velocity- 
dependent  control,  are  also  examined.  Chapter  five  lists  some 
basic  conclusions  about  solar  sail  escape  trajectories.  The 
computer  program,  as  described  in  Appendix  D,  was  the  main  tool 
for  the  analysis  of  solar  sail  escape  trajectories. 


p 


i *' 


2,  FORCE  MODEL 

2,1  General  Force  Model 

The  Intensity  of  sunlight  can  be  expressed  matheiaatlcally 
as  follow3(ref  6)j 

P<  = j ^iT 


where 


irtiere 


P - solar  Intensity  (watt/m^) 
h - Plank’s  constant  (joule-sec) 
iT  - frequency  of  a photon  (1/sec) 

number  of  photons  per  unit  area  per  imlt  time  per 
unit  frequency  (1/m  ) 

Experimentally  this  can  be  expressed  as(ref  5): 


^ - /3S3.^ 


R^-  mean  radius  of  Earth's  orblt=l  a.u.  =1, 495*10 ^^m 
H - distance  to  the  Sun  (m) 

The  solar  pressure,  p,  can  then  be  determined  ast 


p 

-c ' 

^ _ 

(1) 

L • -I 

where 

p - solar  pressure ( newt on/m^) 
c - velocity  of  light  = 3*10®  m/sec 
The  characteristics  of  a reflecting  surface  may  be  defined 
in  terms  of  the  following  parameters: 

r - reflectivity  of  surface 
a - absorption  of  surface 
'2'  - transmissivity  of  surface 
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8 - speculeur  reflection  coefficient 


* emmlsslvlty  of  the  front  and  back  surface 
R g - non-Lambert Ian  coefficients,  front  and  back  surface 

j 

where  Since '^=0  for  the  solar  sall*s  reflecting 

surface,  a=l-r. 

When  light  strikes  a surface,  three  types  of  forces  are 
created  due  to  absorption,  reflection  euid  re-radlatlon  of 
photons.  The  total  force  can  be  expressed  as(ref  6)t 


r = F v-F  (-  f 


where 


F - total  force  (newtons) 

P - force  due  to  reflected  photons  (newtons) 

— \ 

P - force  due  to  absorbed  photons  (newtons) 

“0. 

P - force  due  to  re-radlated  photons  (newtons) 

By  defining  the  unit  vectors  , n,  and  m.  as  In  flgrire  1, 

the  following  relations  are  obtained: 

•/' 

^ -f-  jovn 

/m  - - ^ Of  ^ = X "■  j22^ 

■P-  ^ P-i^ 

where  A Is  the  area  of  the  sail  and  p Is  the  solar  pressure. 

Since  It  may  be  assumed  that  all  Incident  photons  are  Initially 
absorbed  by  the  reflector,  P can  be  expressed  as 


4^-  ■'>1  ^ /9^o(^ 


The  force  due  to  reflected  photons  is  composed  of  two  parts ^ 
due  to  specular (P^)  and  non-specular (F^J  reflection.  The  speculaly 
reflected  light  causes  a force  in  the  -m  direction.  This  Is 
expressed  as 

The  non-speculaly  reflected  light  can  be  thought  of  as  uniformly 
distributed  over  a hemisphere  on  the  reflected  side (see  figure  2). 
The  net  flux  will  be  In  the  -n  direction,  thereby  causing  a 
force  In  the  direction  of  n.  The  total  non-speculaly  reflected 
hemispheric  power  Is  given  by  c(l-s)r(o*A) . Assuming  that 
Lambert’s  Law  applies  and  using  B^to  Indicate  that  the  sail’s 
surface  Is  not  an  Ideal  Lambertian  surface,  the  force  due  to 
non-speculaly  reflected  light  can  be  expressed  as(ref  6)t 

Combining  the  above  expressions,  the  following  equation  for  P 
can  be  obtained) 

- rs 

9 ' 


A 


I 

i ' The  power  emitted  from  a tinlt  area  of  a surface  Is  given 

! by  6<»-T^  where  Is  the  emmlsslvlty  of  the  surface,  Is  the 

I 

Stefan-Boltzman  constant,  and  T Is  the  temperature  of  the  surface. 
Lambert’s  Law  states  that  the  force  per  unit  area  due  to  thermal 

i 

re-radlatlon  on  an  Ideal  Lambertian  siurface  Is  where  n 

Is  the  outward  normal  of  the  surface.  Using  B,. and  B„to  Indicate 

4-  6 

that  the  solar  sail  Is  not  an  Ideal  Lambertian  surface,  the 
following  expression  for  P^can  be  obtained: 

••  — r ' ^ 

Figure  3 shows  a small  area  of  the  sail  that  Is  In  a steady 
state  condition.  Since  the  power  In  has  to  equal  the  power  out, 
the  following  expression  can  be  obtained(ref  6): 
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Dividing  each  side  of  the  equation  by  c and  rearranging  terms: 

Subsltutlon  of  this  result  Into  the  expression  for  F gives: 

“r 

£ -../}(, (5) 

Equation  5 gives  the  Impression  that  It  would  be  best  to 

build  a solar  sail  with  a low  6^.  Upon  Inspection  of  equation 

4,  one  realizes  that  a Is  the  only  parameter  a designer  has  for 

9 

regulation  of  the  temperatiire  of  the  sail,  since  the  character- 
istics of  the  front  side  will  probably  be  fixed  such  that  maximum 
reflection  occurs. 

Combining  equations  2, 3, and  5 gives  the  following  general 
force  model: 


i = 


^ / /■  t S^(^l  ~ ^ r Ciroai  ^ j 

t-  |j^/-rs3  C^OOi 


I 


For  a perfect  reflector  r=l,  s=l,  and  a=»0;  therefore 
eqviatlon  6 reduces  to 

f = CL  (7) 

Equation  6 was  derived  assuming  a flat  sail  with  area  A. 

A similar  differential  force  model  could  be  obtained.  It  must 
also  be  remembered,  with  reference  to  equation  6,  that  all  of 
the  sail's  surface  characteristics  are  functions  of  the  temper- 
ature of  the  sail,  the  frequency  of  the  radiation,  and  the  sur- 
face material  used.  If  the  operating  point  Is  fairly  constant, 
these  values  can  be  approximated  by  constants. 

2.2  JPL  Force  Model 

The  Jet  Propulsion  Laboratory,  using  a differential  form 
of  equation  6,  Integrated  the  force  over  the  surface  of  a sail 
which  may  be  used  for  future  missions.  Angles  9 and  0 were 
defined  as  shown  In  figure  4.  The  Incident  angle  c<  was  veurled 
and  a curve  fit  was  performed  on  the  data.  The  force  on  the 
sail  can  be  expressed  as: 


where 

F - force  f miction  magnitude  (newtons) 

^ - vmit  vector  In  the  direction  of  thrust  vector 
The  force  ftinctlon  magnitude  can  be  expressed  as(ref  l6): 

£ =-  H C^a»<iZQ  ^ C.J  (8) 

where 

^o  *PA 

© - cone  angle  of  the  thrust  vector  with  respect  to 
the  Sun-spacecraft  line. 
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c^jC^fCj  - sail  force  equation  constants 
Prom  figure  4 the  following  expression  can  be  obtained: 

0 

where 

- cenerllne  angle, l.e.  the  angle  between  the  surface 
normal  and  the  thrust  vector,  measured  In  the 
direction  from  n to  1 * . 
e<  - Incident  angle  as  defined  In  figure  1 
A relationship  exists  between  S and  <;ii  , which  Is  shown 
graphically  In  Appendix  A,  along  with  a plot  of  equation  8, 

If  a certain  s Is  desired,  the  relationship  between  <9  and 
Is  used  to  obtain  , thereby  determining  o<  . Knowing  o<  , the 
direction  of  the  normal  of  the  sail  Is  fixed.  The  vectors  1* , 

P,  and  n are  all  In  the  same  plane. 


A perfect  reflecting  sail  would  be  represented  as: 


') 


(9) 


This  Is  equivalent  to  equation  7 


Equation  8 will  be  called  the  JPL  reflector  and  equation  9 
will  bo  called  the  perfect  reflector. 

The  direction  of  the  force  function  can  be  represented  ass 


I 

^ CoO  © it.  ^ 0 Aiyrt  ^ ^ i-  Q 


where  the  cone  8mgle(e)  and  clock  angle (f)  are  defined  In 
figure  5. 


In  order  to  compare  the  performance  of  different  sails,  a 
term  , the  characteristic  acceleration,  Is  defined  as 
follows  s 


Therefore  the  acceleration  of  the  sail  Is 


dr  = 7*  C, 


where 


o 

- acceleration  (m/sec  ) 


^ - 

ao  ■ 


total  mass  of  solar  sail  vehicle  (kg) 

6 


/W1. 


^.Sj3  x/C> 


/ ^cY 

y 


Therefore 


~ <2^  (f/i-Cj  coosB  *■ 


Cj  ^0 


’V 


(U) 


The  maximum  variation  of  a^  due  to  variation  of  R(see 
definition  above)  at  the  sphere  of  Influence  of  Earth  Is  approx- 
imately \%*  Since  this  paper  Is  concerned  with  Earth  escape 
maneuvers  that  are  well  within  the  sphere  of  Influence  of 
Earth,  It  Is  reasonable  to  assume  a^  constant  over  the  maneuver. 
2.3  Effects  of  Solar  Wind  and  Earth  Shine 

The  only  force  acting  on  the  solar  sail  has  been  ass\med 
to  be  the  force  due  to  solar  radiation  pressure.  Other  forces 
that  might  be  considered  are  the  forces  from  solar  wind  and 
Earth  shine. 

The  solar  wind  Intensity  has  been  experimentally  determined 
as  follows (ref  5)* 


The  Intensity  of  can  Increase  by  two  or  three  orders  of  mag- 
nitude during  periods  of  solar  activity.  Comparing  and  ^ 

It  Is  evident  that  solar  wind  pressure  Is  at  least  three  to  four 
orders  of  magnitude  less  than  photon  prossxjre.  Solar  wind  will 
be  neglected  In  this  thesis.  The  major  effect  of  solar  wind  Is 
to  deteriorate  the  sails  surface  and  thereby  decrease  the 
reflectivity. 
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To  determine  the  effects  of  Earth  shine  an  argument  similar 
to  the  one  developed  In  Vlllers  (ref  15)  Is  presented.  The  ab- 
sorption of  sunlight  energy  by  the  Earth  Is  given  as: 

--  Ps-rrRli,.r,) 

where 

E^  - power  absorbed  by  the  Earth  (watts) 

Hg  - radius  of  the  Earth  (6378.l6*lo3  m) 

r^  - albedo/ref lectlvlty  of  the  Earth  (.32  to  .52) 

- solar  Intensity  (watt/m^) 

Assuming  that  all  Earth  shine  effects  occur  at  the  Earth* s 
surface  and  that  the  Earth  re-radlates  Isotropically: 

where 

E^f-  sunlight  power  re-radlated  I'rom  Earth  (watt) 
Earth  re-radlated  Intensity  at  the  Earth’s 
s\irface  (watt/m^) 
substitution  gives 


The  maximum  possible  Earth  shine  Intensity (P^^)  Is  then 


4 

m m 

ti 

5’(-^ 

At  a distance  D above  the  Earth’s  surface  the  energy  recleved 
from  Earth  shine  has  to  be  less  than  or  equal  to  the  maximum 
power  emmltted  at  the  Earth’s  surface. 
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This  paper  assumes  an  operating  point  of  at  least  10  Earth 
radii.  The  effects  of  will  be  at  least  two  orders  of  magni- 
tude less  than  and  therefore  will  be  neglected. 
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3.  OPTIMIZATION 

3»1  Coordinate  Systems 

An  Inertial  coordinate  systemCi)  and  a rotating  coordinate 
system(H)  are  defined  as  in  figure  6 with  the  following  assumptions 
and  definitions* 

1.  Launch  occurs  at  t=0. 

2.  The  inertial  and  rotating  coordinate  system  are 
centered  at  Earth  with  the  x,  y axes  in  the  ecliptic 
plane  and  the  z axis  perpendicular  to  the  ecliptic. 

3.  At  taO,  the  X axis  of  the  inertial  coordinate  frame 
points  at  the  Sun  and  remains  fixed  with  respect  to 
inertial  space  throughout  the  mission. 

4.  The  rotating  coordinate  system  has  its  x axis  pointing 
at  the  Sun  through  out  the  mission. 

5.  Earth's  orbit  about  the  Sun  is  circular. 

6.  The  difference  between  the  Sun-spacecraft  distance 
and  the  Sun-Earth  distance  is  negligible. 
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The  angle  y as  defined  in  figure  6 can  be  expressed  as 


where 

Y - angle  of  Sun  with  respect  to  the  Inertial 
coordinate  frame  (rad) 

^ - rotation  rate  of  Earth  about  the  Svin  ( 1.991 3859*10 
^ sec 

t - time  (sec) 

A rotation  matrix  that  transforms  the  mathematical  repre- 
sentation of  a vector  In  the  R frame  to  Its  equivalent  vector 
representation  In  the  I frame  can  be  defined  as: 
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3.2  Constraints 

The  accelerations  on  the  soleu:  sail  are  caused  by  gravity 
and  solar  pressure.  The  differential  equations  of  motion  can 
be  written  as: 


(15a) 

(15b) 


where 


^ - Earth’s  gravitational  c ons tan t( 398601,2  km3/sec^) 

^ CO 

V - velocity  vector  (km/seo) 
radius  vector  (km) 

r a,  (i-’A’j 
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The  solar  sail  has  as  control  variables  the  cone  and  clock 
angles (eandy)  of  the  vector  The  clock  angle  Is  unconstrained. 

The  cone  angle  Is  constrained  to  be  between  and  180  degrees. 

The  upper  constraint  on  9 results  from  the  mathematical  defini- 
tion of  B and  . The  lower  constraint  results  from  the  fact 
that  the  force  vector  can  never  point  towards  the  Svm.  The  quantity 
Is  the  limiting  cone  angle.  For  the  perfect  reflector.  It 
Is  90  degrees  and  for  the  JPL  reflector  It  Is  the  solution  of  the 
following  equation: 


Cl  i-  i-  ^ o 

this  reduces  to  the  following: 
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Therefore  the  constraints  on  the  control  variables  are: 


For  the  JPL  sqviare  sail  as  given  In  Appendix  A,  ^ Is  eqxial  to 
118.851  degrees. 

3.3  The  Optimal  Control  Law 

Since  the  solsnr  sail  extracts  energy  from  a limitless  source 
for  propulsion  and  attitude  control,  a performance  Index  for 
trajectory  optimization  should  not  Involve  the  energy  source. 

There  Is  an  Inherent  assumption  that  the  control  law  developed 
does  not  require  changes  In  the  control  veurlables  In  a period  of 
time  that  cannot  be  provided  for  by  the  attitude  control  system 
of  the  solar  sail.  The  performance  Index  J Is  given  by: 


The  final  objective  of  the  maneuver  is  to  achieve  escape* 


Therefore  the  final  energy  (E  ) has  to  equal  zero. 
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It  is  also  assumed  that  the  state  (r  and  V)  is  given  at 
t=0.  This  can  be  written  asi 


r‘^Co)  = 


(20a) 

(20b) 


The  problem  is  then  to  minimize  J subject  to  the  boundary 
conditions  in  equations  19  and  20  and  the  constraints  in  eqiiations 
15  and  17*  The  method  used  to  achieve  this  optimization  is  to 
define  a Hamiltonian  which  must  be  maximized  in  order  to 
minimize  J, 

The  Hamiltonian  is  defined  ast 


where 


H - Hamiltonian 
£ (km/sec2) 

pf'"*-  position  co-state  vector  with  respect  to  inertial 
space  (sec/km) 

- velocity  co-state  vector  (also  called  primer  vector) 

— y 

with  respect  to  inertial  space(  sec^/km) 


Applying  the  Euler-Lagrange  equations* 
therefore  ^ 
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therefore 


combining  and  rearranging  terms  gives* 
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(22a) 

(22b) 


Equations  22a  and  22b  are  the  co-state  differential  equations, 

whereas  equations  15a  and  15b  are  the  state  differential  equations 

The  force  vector  and  primer  vector  directions  can  be 
. • 

represented  as  (see  figure  7)* 
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(23a) 

(23b) 


since  the  control  variable  ^ has  iio  constraints,  the 
Hamiltonian  can  be  maximized  with  respect  to  f by  setting  the 
partial  of  the  Hamiltonian  with  respect  to  f equal  to  zero. 
The  Hamiltonian  of  concern  is: 
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f - cr-  (24) 

Prom  equation  24  it  is  evident  that  the  force  vector,  the 
primer  vector,  and  the  x axis  of  the  rotating  coordinate  frame 
will  all  be  in  the  same  plane.  Therefore  the  force  vector's 
direction  can  be  expressed  as  (ref  12) t 
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Because  the  control  variable  e Is  constrained  as  expressed 
In  eqxiatlon  1?,  the  Maximum  Principle  must  be  applied  In  the 
determination  of  the  control  law  for  e«  Substitution  of  equation 
24  Into  the  expression  for  the  Hamiltonian  and  reaurranglng  terms 
gives  the  following  Hamiltonian  of  concern: 

H - 

For  the  Ideal  reflector: 

ML 

Q 7^  3 S — yS  — 

therefore 

3 J ^ *■  SS'.^-n 

y 

The  ambiguity  of  the  sign  In  the  above  equation  Is  removed  by 
ixse  of  the  necessary  condition  that  the  second  partial  of  H* 
with  respect  to  Q has  to  be  less  than  zero  at  a point  that  max- 
Imlzes  H' . Substitution  of  the  above  value  of  8 Into 
yields: 

— ^ 9 ceyo ^/3 i-  S' /Mtf  ^ O 
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since  G Is  In  the  second  quadrant,  the  sign  of  the  radical  has 
to  be  negative.  This  same  result  would  have  been  achieved  If 
physical  reasoning  were  used  to  remove  the  ambiguity  of  the 
above  sign.  If  the  sign  was  positive,  %n  & would  be  greater  than 
zero  and  therefore  e would  violate  the  constraint  eqmtlon. 
Therefore  the  cone  angle  control  law  for  a perfect  reflector  Is: 
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Proceeding  In  a similar  manner  for  the  JPL  reflector.  It 
Is  found  that  the  equations  do  not  solve  analytically  In  closed 
form.  A second  order  Newton  Iteration  method  can  be  used  to 
solve  for  the  optimal  B*  The  Hamiltonian  of  concern  can  bo 
written  as: 

o(e)El_e,ff) 

Where 

0(^)  ~ ^ ^ 

E(e,S)-  <^9o(s-/S) 


defining  the  following: 
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Therefore  the  optimal  0 Is  given  by: 


^ e' 
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Is  given  by  equation  26  (27b) 

^ ^ ^/S'O  , j^g  given  in  equation  l6)  (27c) 

Prom  equation  2ll-  and  26,  It  Is  evident  that  expressions 
for  ^ and  o- must  be  obtained.  Prom  equation  23b  the  following 
expression  for  0 can  be  obtained* 


therefore 


^ c^o 


IS)  M)  ' 

where  and  are  the  x and  y components  of  the  primer  vector 
with  respect  to  the  Inertial  coordinate  frame. 

In  a similar  manner  It  can  be  shown  that* 


Prom  equations  26  and  27  It  Is  clear  that yg  must  be  deter- 
mined before  the  optimal  thrust  cone  angle  (0)  can  be  determined 
Since  In  this  thesis  the  optimal  & will  be  determined  and  then 
checked  to  see  If  It  Is  within  the  constraints,  considerable 

work  could  be  avoided  If  the  constraints  In  equation  17  could 
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be  expressed  In  terms  of  /3 , The  optimal  S Is  obtained  when  the 
partial  of  the  Hamiltonian  with  respect  to  S (HP  as  defined 
above)  Is  set  equal  to  zero.  If  It  Is  assumed  that  the  situation 
Is  such  that  (and  therefore = ) Is  the  optimal  0 f the 

following  would  result: 
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since 
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concluded  that: 


/?{ Bi)  - C,  V / Cj  ~ O 

Is  not  equal  to  zero  for  all  possible  , it  can  be 


Therefore 


/?  = SL- 


Upon  Inspection  of  equation  26  It  Is  realized  that  as  ^ 
approaches  180  degrees,  ^approaches  130  degrees.  Therefore  the 
constraints  Imposed  upon  e In  equation  1?  can  be  replaced  by  the 
following: 
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3.4  Boxmdary  Conditions 

The  boundary  conditions  are  now  addressed.  This  thesis 
assumes  that  the  state  Is  completely  described  at  launch  as 
given  In  equation  20.  Boundary  conditions  on  the  co-state  must 
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now  be  determined  based  upon  equation  19. 

From  Appendix  B the  following  conditions  on  the  co-state 
must  exist  at  the  final  time: 
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For  the  solar  sail  escape  maneuver,  the  following  terms  can  be 
identified: 

3=0 
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Substitution  of  these  terms  Into  eqixatlons  31a  and  31b  and 


\ 


t 


eliminating  V from  the  equations  resxilts  In  the  following 
boundary  conditions  for  the  co-states 

r u)  n 


where 


(32a) 

(32b) 
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The  final  time  (t^)  can  be  obtained  from  equation  31o(equatlon 
19). 

3.6  Summary 

The  analysis  of  an  optimal  Earth  escape  maneuver  to  zero 
energy  results  In  a two-point  boundary  value  problem  of  a 12^^ 
order  system.  The  Implementation  of  the  analysis  In  this  thesis 
Is  to  guess  an  Initial  co-state  and  then  by  use  of  a Runge-Kutta 
algorithm,  Integrate  equations  15  and  22  until  the  energy  equals 
zero.  At  each  point  of  the  Integration  the  control  law  Is 
determined  by  equations  24  and  27  In  conjunction  with  equations 
28  and  29.  The  final  co-state  Is  checked  to  see  If  It  satisfies 
equation  32.  If  It  does  not  (to  a desired  accuracy),  an  Iteration 
scheme  Is  developed  that  Improves  the  Initial  guess  of  the 

co-state.  This  Iteration  la  continued  until  the  desired  accuracy 
Is  achieved. 
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A Newton-Raphson  Iteration  scheme  was  used  to  improve  the 
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Initial  guess  of  the  co-state.  A development  of  the  Newton- 
Haphson  Iterator  Is  presented  In  Appendix  C.  The  computer 
program  developed  to  Implement  the  above  analysis  Is  presented 
In  Appendix  D. 


( 
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4.  DISCUSSION 


4,1  Introduction 

The  purpose  of  this  thesis  Is  to  develop  thfe  optimal  control 
law  necessary  for  a ailnlmura  time  escape  from  a high  Earth  orbit 
with  a solar  sail.  The  end  result  Is  a computer  program 
(see  Appendix  D)  that  determines  the  time  history  of  the  optimal 
control  law.  With  minor  changes,  this  computer  program  could 
be  used  to  determine  optimal  escape  trajectories  from  planets 
other  than  Earth, 

The  computer  program  developed  does  not  have  active  con- 
straints to  prevent  the  optimal  trajectory  from  going  below 
a certain  desired  altitude.  This  was  not  a serious  problem  for 
escape  trajectories  whose  Inltal  orbit  was  above  100000  kilometers. 
Some  escape  trajectories  from  Initial  ecliptic  orbits  with  a 
semi-major  axis  of  100000  kilometers  or  less  did  Intersect  the 
surface  of  the  Earth,  A recommended  solution  to  this  problem 
Is  to  add  a penalty  ftinctlon  to  the  performance  Index  which 
heavily  penalizes  trajectories  close  to  the  surface  of  the  Earth, 

A Newton-aaphson  Iteration  scheme  was  used  to  solve  the 
two-point  boundary  value  problem.  The  computer  program  devel- 
oped had  problems  converging  for  some  Initial  orbits.  The  con- 
vergence problems  were  usually  caused  because  the  initial  guess 
of  the  co-state  was  too  far  from  the  optimal  initial  co-state 
for  the  Newton-Haphson  Iterator  to  work.  One  possible  solution 
to  this  problem  Is  to  develop  a gradient  method  for  solving  the 
two-point  boundary  value  problem  with  a velocity-dependent 


control  law  (discussed  below)  as  the  Initial  guess.  Another 
possible  solution  to  this  problem  is  to  use  the  velocity- 
dependent  control  law  to  obtain  an  estimate  of  the  optimal 
escape  trajectory.  Using  equation  J2  and  the  estimate  of  the 
optimal  escape  trajectory,  an  approximation  to  the  final  optimal 
co-state  can  then  be  obtained.  The  computer  program  in  Appen- 
dix D can  then  be  used  to  Integrate  the  state  and  co-state 
backwards  In  time  until  the  initial  or  near  initial  state  is 
reached.  The  initial  co-state  as  given  by  the  backward  integra- 
tion will  usually  be  sufficiently  close  to  the  optimal  initial 
co-state  such  that  a converged  escape  trajectory  can  be  obtained. 
This  latter  method  has  proved  successful  in  this  thesis. 

Although  this  thesis  is  concerned  mainly  with  the  optimal 
control  law,  two  other  sub-optimal  control  laws  were  examined. 
They  are  called  open-loop  control  and  velocity-dependent  control. 
An  open-loop  control  trajectory  results  when  equations  15 
and  22  are  integrated  forward  in  time  with  the  initial  co-state 
vectors  as  given  below: 
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(33b) 

The  quantities  r'*^  and  are  the  initial  state  vectors  with 

respect  to  the  inltal  coordinate  frame,  and  V la  the  circular 
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orbital  velocity  associated  with  r'^.  At  each  point  of  the  inte- 
gration the  control  variables  & and  f are  determined  by  equations 
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24  and  27  in  conjunction  with  equations  28  and  29.  No  attempt 
is  made  In  open-loop  control  to  satisfy  equation  32.  An  open- 
loop  control  trajectory  would  be  the  optimal  control  trajectory 
If  the  Initial  guess  of  the  co-state  (eqviatlon  33)  was  the 
optimal  Initial  co-state.  Equation  33  results  from  the  fact 
that  tangential  thrust  has  been  shown  to  be  a good  approximation 
to  the  optimal  thrust  for  Initial  circular  orbits  (ref  4) . 

Velocity-dependent  control  results  when  the  acceleration 
vector  (equation  11)  Is  required  to  lie  In  the  plane  of  the  r 
axis  of  the  rotating  coordinate  frame  euid  the  velocity  vector 
and  Is  oriented  In  the  direction  that  maximizes  Its  projection 
on  the  velocity  vector.  A velocity-dependent  control  trajectory 
results  when  eqviatlon  15  Is  Integrated  forward  In  time  with  f 
and  0 being  determined  at  each  point  of  the  Integration  by  equa- 
tions 24  and  2?  In  conjunction  with  .7-  and  ^ as  defined  below: 
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Satisfaction  of  the  constraints  In  equation  1?  Is  also  required. . 
Velocity-dependent  control  relies  only  upon  the  state,  and  there- 
fore there  Is  no  need  for  the  co-state  equations  and  boundary 
conditions. 

4.2  Optimal  Control 

Figures  8 and  9 are  the  time  histories  of  the  optimal 


control  variables  for  a minimum  time  escape  from  an  Initial 
150000  kilometer  circular  ecliptic  Earth  orbit.  The  maximum 
rate  of  change  of  the  control  variable  9 la  approximately  one 
degree  per  hour.  The  JPL  square  sail  can  sustain  an  attitude 
rate  of  between  20  to  25  degrees  per  hour.  The  Infinite  atti- 
tude rates  of  the  control  variable  f results  from  the  mathemat- 
ical singularities  associated  with  cone  and  clock  angle  repre- 
sentation. At  1,7  and  4,5  days  9 Is  at  180  degrees  and  there- 
fore !/>  can  have  any  value  for  the  same  force  vector.  Because 
the  JPL  square  sail  Is  not  Ideal,  Its  , as  given  by  equation 
16,  Is  not  90  degrees.  This  Is  Indicated  as  a straight  line 
between  5.9  and  9*1  days  on  figure  8.  During  this  period  y can 
be  rotated  through  the  necessary  180  degrees  while  Q Is  maintained 
at  This  results  In  an  attitude  rate  of  about  2,5  degrees 

per  hour.  Similar  results  were  obtained  for  Initial  100000 
kilometer  circular  ecliptic  orbits  with  resulting  attitude  rates 
about  twice  as  large.  The  attitude  rates  of  the  solar  sail 
tend  to  be  Inversely  proportional  to  the  semi-major  axis  of  the 
Initial  orbit. 

Figure  10  demonstrates  some  of  the  basic  characteristics 
of  optimal  escape  trajectories  from  ecliptic  orbits  with  a solar 
sail.  The  trajectory  tends  to  become  very  eccentric  quickly 
with  the  Instantaneous  semi-major  axis  of  the  trajectory  perpen- 
dicular to  the  Sun-spacecraft  line  In  the  orbital  plane.  The 
Instantaneous  perlcenter  of  the  trajectory  Is  on  the  side  of  the 
Earth  on  which  the  solar  sail  can  achieve  Its  maximum  thrust  In 
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the  direction  of  the  velocity  vector.  Figure  11  shows  the 
Increase  In  eccentricity  for  various  initial  ecliptic  orbits. 

The  optiT^l  control  law  for  an  Initial  orbit  Inclined  to 
the  ecliptic  plane  Is  shown  In  figure  12.  The  Initial  orbit 
considered  Is  the  initial  circular  ecliptic  orbit  in  figure  10, 
rotated  out  of  the  ecliptic  plane  by  JO  degrees.  Symmetry  about 
the  x-z  plane  Is  maintained.  The  control  law  for  e is  similar 
to  figure  8,  The  time  to  escape  Is  0,1  days  longer  than  that 
from  a similar  Initial  ecliptic  orbit.  The  control  law  for  'j' 
is  less  demanding,  and  the  trajectory  tends  to  become  eccentric 
slower  than  It  does  for  a similar  Initial  ecliptic  orbit.  The 
Instantaneous  semi-major  axis  of  the  orbit  is  in  the  orbital 
plane  perpendicular  to  the  Sun-spacecraft  line.  As  the  Inclina- 
tion of  the  Initial  orbit  is  increased,  the  control  law  tends 
to  be  less  dema.ndlng,  and  the  escape  time  tends  to  be  longer. 

The  optimal  control  law  for  an  Initial  150000  kilometer 
circular  polar  orbit  Is  shown  In  figures  13  and  1^.  The  control 
variable  & has  a very  low  attitude  rate.  The  attitude  rate  of 
■f'  Is  about  2.2  degrees  per  hour  at  first  and  then  falls  off 
considerably.  The  optimal  escape  trajectory  from  circular  polar 
orbits  maintains  a low  eccentricity  throughout  most  of  the  escape 
maneuver  as  shown  In  figure  15 . Similar  results  were  obtained 
for  lower  Initial  polar  orbits.  As  can  be  seen,  the  attitude 
rates  for  Initial  polar  orbits  are  less  than  those  for  similar 
Initial  ecliptic  orbits  but  the  final  time  Is  longer. 

It  appears  that  the  performance  of  a solar  sail  will  tend 
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to  be  between  the  limiting  cases  of  the  resulting  escape  trajec- 
tories from  similar  ecliptic  and  polar  initial  orbits.  The 
ecliptic  initial  orbit  will  yield  the  highest  attitude  rates 
that  can  be  expected  and  the  smallest  final  escape  time.  The 
polar  initial  orbit  will  yield  the  lowest  attitude  rates  and 
the  largest  final  time, 

4.3  Open-Loop  Control 

Open-loop  control  is,  as  the  name  implies,  an  open  loop 
trajectory  of  the  optimal  control  law  with  the  initial  guess  of 
the  co-state  vectors  as  given  in  equation  33 • For  initial 
circular  ecliptic  orbits,  open-loop  control  is  a good  approxi- 
mation to  the  optimal  control. 

For  the  initial  I50OOO  kilometer  circular  ecliptic  orbit 
in  flgiire  8,  the  final  escape  time  predicted  by  open-loop  control 
varies  from  the  optimal  control  final  escape  time  by  less  than 
1%,  For  the  initial  100000  kilometer  circular  ecliptic  orbit 
in  figure  11,  the  error  is  around  As  the  eccentricity  of 

the  initial  orbit  is  Increased  the  error  between  open-loop  and 
optimal  control  Increases.  For  the  initial  eccentric  ecliptic 
orbit  in  figure  11  (100000  kilometer  perlcenter,  0.25  eccentricity), 
the  error  in  the  final  escape  time  between  these  two  methods  is 
about  8.5/^*  As  the  inclination  of  the  initial  orbit  with 
respect  to  the  ecliptic  is  Increased,  the  error  between  optimal 
and  open-loop  control  also  increases.  It  can  generally  be 
concluded  that  the  open-loop  control  law  is  a good  approximation 
to  the  optimal  control  law  for  a minimum  time  escape  from  high 


36 


Earth  orbits  with  low  inclination  and  low  eccentricity. 

Figures  l6a  and  l6b  were  constructed  using  open-loop  control. 
Figure  l5a  shows  the  final  escape  time  for  four  various  starting 
positions  fron  a 15OOOO  km  circular  ecliptic  orbit.  Figure  l6b 
shows  how  the  energy  Increases  during  the  escape  maneuver  of  two 
of  the  trajectories  considered  in  figure  l6a.  In  figure  l6b  the 
sharp  Increases  in  energy  followed  by  little  or  no  increases  is 
because  in  the  ecliptic  plane  the  solar  sail  Increases  its  energy 
from  solar  radiation  pressure  only  when  it  is  traveling  away 
from  the  Sun.  From  figure  l6b  it  can  be  seen  that  the  initial  en- 
ergy of  the  solar  sail  determines  how  much  energy  must  be  added  in 
order  to  achieve  escape.  Given  the  initial  energy,  certain  start- 
ing positions  make  it  possible  for  the  solar  sail  to  achieve  es- 
cape prior  to  making  another  revolution  over  that  portion  of  the 
trajectory  where  the  vehicle  is  traveling  towards  the  Sun.  This 
is  the  case  for  ■Y'  equal  to  270  degrees  in  figure  l6b.  For  equal 
to  90  degrees  in  figure  l6b,  the  effects  of  having  to  traverse  the 
•backside**  of  the  trajectory  an  additional  time  in  order  to  reach 
escape  condition  is  obvious.  It  can  also  be  seen  that  if  the 
eccentricity  of  the  orbit  is  changed,  the  initial  energy  and  start 
Ing  position  remaining  the  same,  the  final  escape  time  will  bs 
affected  because  the  time  in  which  the  solar  sail  is  in  "useful" 
sunlight  is  changed.  It  therefore  appears  that  the  final  escape 
time  from  ecliptic  orbits  is  dependent  not  only  upon  the  initial 
energy  of  the  vehicle  but  also  upon  its  starting  position  within 
the  orbit  and  the  eccentricity  of  the  initial  orbit, 
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4.4  Veloclty«Dependent  Control 

The  velocity-dependent  control  law  orients  the  force  vector 
of  the  sail  In  a direction  that  mazlmlzes  Its  projection  on  the 
velocity  vector  In  the  plane  of  the  velocity  vector  and  the  Siin- 
spacecraft  line  while  still  satisfying  equation  17.  Velocity- 
dependent  control  Is  rather  interesting  In  that  It  Is  slapller  to 
Implement  and  analyze  and  yet  gives  results  that  are  close  to 
optimal. 

For  the  Initial  150000  1cm  circular  ecliptic  orbit  In  figure 
8,  velocity-dependent  control  yields  a final  time  of  20,33  days. 

The  optimal  final  escape  time  was  20,29  days,  and  the  open-loop 
control  law  yielded  a final  time  of  20,29  days.  This  gives  an 
error  of  about  1.2^^  between  velocity-dependent  and  optimal  control. 
For  the  100000  km  circular  ecliptic  Initial  orbit  In  figure  11, 
the  error  In  final  time  between  velocity-dependent  and  optimal 
control  was  1,?^,  As  the  Initial  orbit  Is  Inclined  to  the  ecliptic 
plane,  the  error  between  velocity-dependent  and  optimal  control 
Increases,  For  the  Initial  I50000  km  circular  polar  orbit  In 
figures  13  and  14,  the  optimal  escape  trajectory  gave  an  escape 
time  of  29.45  days.  The  velocity-dependent  control  law  yielded 
a final  escape  time  of  30.42  days.  This  Is  an  error  of  about 
3.3^  between  velocity-dependent  and  optimal  control.  The  open- 
loop  control  law  gave  a final  escape  time  of  32.7  days  which 
results  In  an  error  of  11%  from  the  optimal  control  law. 

The  velocity-dependent  control  law  was  used  to  obtain 
flgurej*17a  and  17b.  Figure  17a  shows  that  starting  position 
In  the  ecliptic  plane  has  little  effect  upon  the  final  escape 
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time.  Figure  1713  shows  that  for  a polar  orbit  the  Increase  In 
energy  Is  nearly  uniform  through  out  the  escape  maneuver.  There- 
fore the  final  escape  time  from  polar  orbits  Is  almost  totally 
dependent  upon  the  Initial  energy  of  the  vehicle.  Starting 
position  In  the  orbit  and  the  eccentricity  of  the  Initial  orbit 
have  little  effect  upon  the  escape  time  from  polar  orbits 

The  velocity-dependent  control  law  appears  to  be  a good 
approximation  for  the  optimal  control  law  of  a minimum  time 
escape  from  high  Earth  orbits.  As  the  orbit  Is  Inclined  the 
velocity-dependent  control  law  Is  a better  approximation  to  the 
optimal  control  law  than  open-loop  control. 


Eooentrlolty  versus  Time  for  various  Initial  orbits  - Optimal  Control 


Optimal  Control  - Theta  versus  Time 


Initial  circular  polar  orbit  150000km  k 

1.63014  km/sec  2. 

JPL  square  sail  ao=  1 mm/sec2  “* 


1 1 

I 

1 

J 

I 

optimal  Control  - Eccentricity  versus  time 
Initial  circular  polar  orbit  - 15000  km  k ' 

: -1.63014  km/sec  J. 

JPL  sq’iare  sail  - a^  = 1 mm/sec^ 


Plgitre  15 
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Final  Time  versus  Different  Starting  Positions 


velocity-dependent  control  JPL  square  sail  a^=  1 mm/sec 

circular  polar  orbits  with  r = 150 CIO U km 

o 

angle  in  the  orbital  plane 
as  shown  below:  . .W 


1 


Figure  17b 
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CONCLUSIONS 


1.  Escape  from  a high  SeLrth  orbit  with  solar  sail  propulsion 
Is  possible. 

2.  Optimal  escape  trajectories  In  the  ecliptic  plane  will  give 
the  minimum  escape  time,  but  will  demand  the  greatest  attitude 
rates.  Ecliptic  escape  trajectories  Involve  very  fast  Increases 
In  eccentricity,  with  semi-major  axes  perpendicular  to  the  Sun- 
spacecraft  line.  Escape  time  from  ecliptic  orbits  Is  dependent 
upon  the  Initial  energy  of  the  spacecraft,  the  starting  position 
of  the  vehicle  In  the  orbit,  and  the  eccentricity  of  the  Initial 
orbit. 

3.  Optimal  polar  escape  trajectories  will  have  the  least  demanding 
attitude  rates,  but  the  escape  times  will  be  longer  than  similar 
Inclined  and  ecliptic  escape  trajectories.  Polar  escape  trajectories 
maintain  a low  eccentricity  throughout  most  of  the  trajectory. 

Escape  time  from  polar  orbits  Is  dependent  upon  the  Initial 

energy  of  the  spacecraft.  Starting  position  In  the  orbit  and 
the  eccentricity  of  the  Initial  orbit  have  little  effect  upon 
the  final  escape  time  for  polar  escape  trajectories. 

4.  The  performance  of  optimal  escape  trajectories  Inclined  with 
respect  to  the  ecliptic  plane  will  tend  to  be  between  the  limit- 
ing cases  of  similar  polar  and  ecliptic  escape  trajectories. 


5.  Open-loop  control  Is  a good  approximation  to  the  optimal 
control  law  of  a minimum  time  escape  with  a soleo*  sail  from 

a high  Earth  v.  )lt  with  low  Inclination  and  low  eccentricity. 

6.  Open-loop  control,  with  the  proper  scale  factor,  serves  as 

a good  initial  guess  for  the  co-state  vectors  In  the  determina- 
tion of  the  optimal  control  law.  The  scale  factor  should  be 
chosen  such  that  the  final  co-state  Is  of  the  same  order  of 
magnitude  as  the  calculated  optimal  final  co-state  (equation  32). 

7.  Velocity-dependent  control  serves  as  a good  approximation 
for  the  optimal  control  law  of  a minimum  time  escape  from  a- 
high  Earth  orbit. 
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APPENDIX  A 


Force  Equation  Parameters 


The  values  of  all  constants  used  In  this  thesis  are: 


Perfect  Sail 

JPL  Square  Sail 

JPL  Helloff-yro 

r 

1 

0.38 

0.83 

s 

1 

0.9^ 

0.94 

0 

0.05 

0.05 

^6 

0 

0.55 

0.55 

2/3 

0.79 

0.79 

2/3 

0.55 

0.55 

0.5 

0.3^9 

0.367 

0.5 

0.662 

0.643 

c 

3 

0 

-0.011 

-0.010 

The  values  for  the  constants  In  the  JPL  square  and  hello- 
gyro  sails  were  obtained  from  a JPL  Interoffice  memorandum  dated 
28  March  1977  (ref  l6).  These  are  the  proposed  specifications 
on  possible  futxire  solar  sails  being  considered  by  the  Jet 
Propulsion  Laboratory.  At  high  altitudes,  where  the  attitude 
controls  are  slowly  varying,  the  hellogyro  sail  can  be  approxi- 
mated by  the  model  developed  In  this  thesis. 

Typical  values  of  a^  for  the  JPL  square  and  hellogyro  sails 
vary  from  about  0,3  to  1,2  mllllmeters/sec^ , This  thesis  assumes 
an  a^  of  1 mm/sec^  for  all  calculations. 
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Pljarure  A-1  Is  a graph  of  the  normalized  force  versus  B for 
the  three  types  of  solar  sails  listed  above.  Equation  8 is  the 
generating  function.  The  quantity  6 Is  the  cone  angle  of  the 
sail's  force  vector  with  respect  to  the  x axis  of  the  rotating 
coordinate  frame  (see  figure  5), 

Figures  A-2  and  A-3  are  polar  plots  of  equation  11  for  the 
ideal  and  JPL  square  sail.  The  polar  plots  are  In  the  plane  of 
the  primer  vector  and  the  Sun-spacecraft  line.  The  only  quadrant 
of  concern  Is  the  upper  right  quadrant  since  this  Is  were  G 
satisfies  equation  17,  The  upper  constraint  on  the  cone  angle 
In  equation  17  (180  degrees)  results  from  the  mathematical  defin- 
ition of  a cone  angle.  The  lower  constraint  In  equation  17  {oj) 
results  from  the  physical  fact  that  the  sail's  force  vector  can 
never  point  toward  the  Sun.  From  figure  A-1  It  can  be  seen  that 
6^  Is  around  118  degrees.  Figures  A-2  and  A-3  also  show  the 
basic  principle  behind  the  optimal  d»  Given  the  direction  of 
a primer  vector,  the  optimal  b Is  the  angle  that  maximizes  the 
projection  of  P onto  When  ^ Is  less  than  , the  maxlmxim 
projection  of  P onto  results  when  F equals  zero  or  & equals 
Figure  A-4  Is  a plot  of  the  centerline  angle  («;»  versus 
the  cone  angle  (0).  The  centerline  angle  Is  the  angle  between 
the  sail's  normal  and  the  thrust  vector  (see  figure  4),  The 
smalysls  developed  In  this  thesis  determines  the  optimal  9 . 

Using  figure  A-4  or  an  analytical  representation,  the  correspond- 
ing optimal  0 can  be  determined.  The  following  equation  can  be 
used  to  find  the  Incident  angle  o<,(see  figure  4): 
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The  analysis  developed  in  this  thesis  also  determines  the  optimal 
V'  of  the  force  vector  (see  flgvire  5).  Knowing  and  , the 
orientation  of  the  solar  sail  Is  defined  and  the  sail  can  then 
be  positioned  In  space  to  achieve  the  optimal  thrust. 


Polar  Plot  of  Normalized  Aocereratlon  Vector  for  JPL  Square  Sail 


polar  plot  In  the  plane  of  the  primer 
vector  and  the  Sun-spacecraft  line 


a, 

normalized  acceleration  vector  = "Si; 
the  arrow  on  the  cvirve  Indicates  the 
direction  of  the  normalized  acceleration 
vector  as  theta  goes  from  90  degrees  to 
180  degrees.  When  theta  Is  less  than 
the  acceleration  vector  Is  negative. 


Figure  A-3 
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APPENDIX  B 


W 


Optimal  Control  Theory 


This  appendix  develops  the  basic  theory  of  optimal  control 
used  In  this  thesis  (ref  2 and  13). 

Let  the  Integral  I be  defined  as: 


2"=  / r(2L,x,±)  Jit 


(B-l) 


If  the  problem  Is  to  find  the  function  x that  maximizes  I,  the 
Integrand  must  satisfy  the  Euler-Lagrange  equation  (ref  13); 


^ _ J.  iL 

Jdt  ^2L 


(B-2) 


Also  If  the  end  point  Is  assumed  to  be  fixed  (x^  and  ti 
specified) , the  following  generalized  boundaury  condition  must 
be  satisfied  (ref  13): 


4F-(ji)aJ  St 


^ o 


(B-3) 


The  problem  that  this  thesis  is  concerned  with  Is  the 
minimization  of  the  performance  Index  Ji 


where 


(bU) 
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S=S(x(tf),tf) 

L=L(x,u,t) 

x=state 

u=control 

tfsflnal  time  (free) 
tj^=lnltlal  time  (fixed) 

It  l3  also  assximed  that  the  folloifing  constraint  equations  are 


imposed: 


(state  equations) 


(functions  specified 
at  final  time  ) 


Using  the  method  of  Lagrange  multipliers,  the  constraint  equations 
can  be  adjoined  to  J,  and  therefore  the  redefined  performance 
index  J*  can  be  treated  in  an  vinconstralned  problem.  Therefore: 

A -r 

J f J [t  i-  i 


where 


-S  t^)-  Si-  !l_ 

- Lagrange  multiplier  for  constraint  at  final  time, 

A ” Lagrange  multiplier  for  state  equations  (co-state). 


The  following  performance  index  can  be  defined: 
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since  x(tj^)  and  tj,  are  assumed  fixed  the  minimization  of  PI  Is 
the  same  as  the  minimization  of  J*.  Therefore  J*  can  be  re- 
vrrltten  as: 


^ * hfJ  ^ 

By  use  of  the  chain  rule  this  can  be  re-wrlttem  as; 

Defining  the  Hamiltonian  H such  that  maximizing  It  minimizes  J; 


i)  = -L 


(B-5) 


Therefore  J*  can  be  re-wrltten  as 


_ ^ _ A 

t. 


(B-6) 


Applying  the  Euler -Lagrange  equation  to  J*  and  remembering  that 
X In  equation  B-2  corresponds  to  x and  u In  equation  B-6 


iil  _ ^1^1  _ ^ 

J±  [_J*J  ^ 

_ ^ ^ ^ A/  - ^ 


f 


Ifl- 


(B-7) 


(B-S) 


from  Inspection  of  the  Hamiltonian,  the  following  expres- 
sion can  be  obtained; 

^ N 


(B-9) 


Applying  the  generalized  boundary  condition  (equation  B-3), 
the  following  expression  can  be  written: 

JL' 


^2. 


*■ 


■ - (-&ji 


Substitution  of  L*  gives  the  following: 


42. 


Pi'\  • ii  U /■>!■'  \T\  . 

(.Jij-  *■  H ~ // '■Aa.  -(jl_ 


Itl 


-o 


Substitution  of  S’  and  H and  rearranging  terms  gives: 


[iS  , .r 


ji? 


'ire. 


therefore 


^ _ ii  il 


nT 


t - 


H Ji 


y ~ 


Ji 


^ Q 


t:t. 


=s 
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(B-lOa) 

(B-lOb) 

(B-lOe) 


Therefore  In  order  to  minimize  equation  a Hamiltonian, 

defined  as  In  equation  B-5»  must  be  maximized.  Equations  B-7 , 
B-3  and  B-9  give  the  differential  equations  and  control  law  that 
must  be  satisfied,  subject  to  the  boundary  conditions  on  the 
state  at  the  initial  time  and  the  co-state  at  the  final  time 
(equation  B-10),  in  order  to  determine  the  optimal  escape 
trajectory. 
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APPENDIX  C 


An  Energy  Dependent  Newton-Haphson  Iterator 


The  solution  to  the  problem  of  minimizing  the  final  time 
on  an  Earth  escape  maneuver  to  zero  energy  has  been  shovm  to  be 
the  solution  of  the  following  set  of  coupled  differential 
equations: 


(C-la) 


(C-lb) 


(C-lc) 


(C-ld) 


Tne  optimal  control  variables,  £)  and  </ , are  determined  by  equations 
24  and  27  In  conjunction  with  equation  28  suid  29,  The  boundary 
conditions  for  the  above  set  of  differential  equations  are: 


L°) 


(C-2a) 


)l%)  - lo 


(C-2b) 


(C-2c) 
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(C-2d) 


and  * denotes  the  optimal  end  conditions  as  developed  in  Appen- 
dix B.  For  an  escape  trajectory  Ef  is  zero. 


As  can  be  seen,  this  is  a two-point  boundary  value  problem. 
The  method  of  solution  will  be  to  guess  and  and  then 


Integrate  equation  C-1  forward  until  equation  C-2e  is  satisfied, 
rne  results  of  equations  C-2c  and  C-2d  are  then  compared  to  the 
final  co-state  of  the  nximerlcal  integration.  If  the  two  co- 
otates  agree  to  within  a desired  accuracy,  the  problem  is  solved. 
If  they  do  not  agree,  an  iteration  scheme  will  be  used  to  Improve 
the  initial  guess  on  the  co-state.  This  iteration  scheme  is 
developed  as  follows. 

Define  the  following  function  h* 

A - & - A 

where 

(from  numerical  integration) 

(from  equation  C-2c  and  C-2d) 
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If  it  Is  assumed  that  as  E(x(t))  approaches  Ef,  energy  Is 
a monotonlc  increasing  function,  a one-to-one  mapping  will  exist 
between  t^,  and  Ef.  In  other  words  If  a given  Ef  is  specified 
for  a given  initial  state  and  co-state,  the  final  time  is 
uniquely  specified.  If  Ef  la  changed,  tf  changes  correspondingly 
It  is  also  evident  that  if  Ef  is  held  fixed  and  the  initial 
co-state  is  varied,  tf ,x(tf ) ,^(tf ) and  h will  all  vary.  The 
quantity  h is  therefore  a function  of  the  initial  co-state  and 
the  specified  terminating  energy  condition  for  a fixed  initial 
state.  This  can  be  represented  as: 

k ^ ll(f(.o),£^)  (0-5) 

By  use  of  a Taylor’s  series  expansion,  the  following  is 
true  to  first  order: 


For  the  problem  defined  in  this  thesis,  Ef  is  fixed;  there- 
fore equals  zero.  Also  the  desired  change  in  should 

be  such  that  to  first  order,  h is  driven  to  zero  (final  co- 
state and  optimal  final  co-state  are  equal).  Therefore: 

Q = k(f(a4)->- 

Define  the  matrix  G as  follows:  ^ 


Hear ranging  terms  and  substitution  gives* 


(C-6) 


The  improved  guess  on  the  initial  co-state  can  be  written  as* 


(C-7) 


The  implementation  of  equation  C-6  requires  the  determin- 
ation of  the  matrix  G.  Analytically  this  is  very  difficult. 

The  method  used  in  this  thesis  is  to  calculate  neighboring  tra- 
jectories. The  development  of  equation  C-5  has  shown  that  this 
problem  can  be  considered  in  an  "energy  space"  Instead  of  the 
conventional  "time  space".  Each  neighboring  trajectory  should 
therefore  be  terminated  when  the  energy  equals  Ef.  The  result 
of  this  technique  is  that  G is  a six-by-six  matrix  and  six 
neighboring  trajectories  need  to  be  calculated  on  each  iteration. 

An  obvious  difficulty  of  working  in  a so-called  "energy 
space"  is  that  equation  C-1  has  to  be  Integrated  in  a "time 
space".  Since  the  Eunge-Kutta  integration  scheme  takes  finite 
time  steps,  satisfaction  of  the  energy  stopping  condition  will 
generally  occur  between  Integration  steps.  A simple  interpola- 
tion is  used  to  adjust  the  size  of  the  last  Integration  step. 

This  is  shown  graphically  in  figure  C-1. 

VThen  the  energy  becomes  greater  than  Ef  on  a particular 
step,  the  average  slope  of  the  energy  curve  over  the  last  step 
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Is  determined  as 


(C-8) 


The  stepslze  ss^,  which,  If  taken  from  tQ,  will  yield  an  energy 
approximately  equal  to  Ef  at  tato+ss^,  can  be  expressed  asx 


ss 


lO 


(C-9) 


The  state  and  co-state  are  re-lnltallzed  to  their  values  at 
t=to,  and  the  time  step  ss^  Is  taken  so  that  condition  E=Ef 
Is  approximately  achieved. 
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APPENDIX  D 


w 


Computer  Program 

SOLCO  (Solar  Control)  Is  the  name  of  the  computer  proaram 
developed  to  determine  the  optimal  minimum  time  escape  trajectory 
of  a spacecraft  propelled  by  solar  radiation  pressure.  SOLCO 
solves  the  two-point  boundary  value  problem  by  use  of  a Newton- 
Raphson  iterator  (see  Appendix  C).  The  sensitivity  matrix  needed 
In  the  Newton-Raphson  iterator  Is  determined  by  comparlns; 
neighboring  trajectories  to  the  nominal  trajectory. 

SOLCO  Is  vrrltten  In  Fortran  IV  and  uses  double  precision 
arithmetic.  SOLCO  requires  a central  memory  core  of  about 
60*10^  bytes  and  a CPU  time  of  between  one  to  six  minutes  on 
an  Amdahl  470  computer.  The  CPU  time  depends  mainly  on  the  Initial 
orbit  of  concern  and  the  guess  for  the  unknown  Initial  co-state. 

SOLCO  consists  of  a MAIN  PROGRAM  and  9 subprograms.  The 
MAIN  PROGRAM  and  each  of  the  subprograms  are  discussed  separately.  ! 

MAIN  PROGRAM 

This  Is  the  controlling  orogram  for  SOLCO, 

I 

1 

! 

Subprogram  INPNCL  j 

This  subprogram  Is  called  by  the  MAIN  PROGRAM,  All  Initial  i 

I 

data  Is  read  and  constants  set  In  this  subprogram.  The  Initial  i 

i 

data  and  constants  are  also  Printed  out.  This  initial  data  Is  I 

I 

read  In  MKS  units  and  converted  to  reference  orbit  canonical  | 

71  i 


units  vhere  1 DU=6373.l6  ka,  1 rU=S06.3l59  sec  and  /<  =1  DU3/TU2. 
The  data  Is  read  in  with  respect  to  the  Inertial  coordinate  fraae 
as  discrlbed  In  Chapter  3*  The  input  variables  are  described 
below, 

AC  = chacterlstJLc  acceleration  (mm/sec^) 

y(l)  - y(3)  = radius* vector  (ka) 

y(^)  - y(6)  = velocity  vector  (km/sec) 

y(7)  - y(9)  = position  co-state  vector  (sec/ka) 

y(10)  - y(l2)  = velocity  co-state  vector  (sec^/km) 

prmtd)  = initial  time  (days),  usTially  set  equal  to  zero. 

prat (2)  = terminating  time  (days)  usually  set  large  enough  so 

that  it  does  not  interfer  with  the  energy  terminating 
condition, 

prnit(3)  = integration  stepsize  (days).  For  an  orbit  with  a 

semi-major  azls  of  100000  km,  a stepsize  of  0.08  day 
gave  good  results.  For  a semi-major  axis  of  I5OOOO  km, 
0.01  day  was  used. 

EF  = terminating  energy  condition  (km^/sec^).  Usually  set  to 
escape  condition,  EF=0. 

IROT  = Earth  rotation  factor.  If  the  effects  of  the  Earth's 
rotation  about  the  Sun  is  desired,  IHOT  should  not 
equal  zero, 

IPH  3 Print  step.  The  value  of  IPR  causes  every  IPH^h  output 
of  the  Runge-Kutta  subprogram  to  be  printed;  that  is, 

IPR=1  causes  every  output  to  be  printed. 

INP  =»  Plotting  factor.  This  program  has  the  option  of  storing 
the  data  from  each  print  step  on  tape/file  so  that  it  can 
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be  used  later,  INP  equal  to  a non-zero  number  will  cause 
storage. 

KSTEP  = Stepsize  factor,  KSTEP  determines  whether  the  nelghbborlng 
trajectories  in  subprogram  NEWIT  are  perturbed  by  a frac- 
tional amount  of  the  nominal  values  (KSTSP=0)  or  by  a 
fixed  amount  (KSTSP=1).  KSTEP=0  was  usually  used, 

NIMAX  =»  Maximum  number  of  iterations  allowed  in  NEWIT, 

FLIM  * desired  accviracy  of  the  Newton-3aphson  iterator.  FLIM  = 

0,1  gave  good  results, 

xs(l)  - xs(6)  =!  stepsize  of  the  neighboring  trajectories  of  the 

co-St  te  vectors.  Usually  taken  as  10“^  with 
KSTEP=0 , 


Subprogram  NEWIT 

This  subprogram  is  called  by  the  MAIN  PHOGRAM,  This  is  the 
Kewton-Haphson  iterator  as  dlscused  in  Appendix  C,  A nominal 
trajectory  is  calculated.  Errors  in  the  final  co-state  conditions 
are  determined.  The  initial  co-state  vectors  are  varied  and 
neighboring  trajectories  are  calculated  in  order  to  get  a six- 
by-six  sensitivity  matrix.  This  matrix  is  inverted  and  post- 
multiplied  by  the  final  co-state  error  in  subprogram  DCHOUT. 

New  values  of  the  initial  co-state  vectors  are  obtained,  and  a 
new  nominal  trajectory  is  calculated.  This  is  continued  until 
the  magnitude  of  the  error  in  the  final  co-state  is  less  than 
some  desired  accuracy  (FLIM)  or  until  convergence  fails. 
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Subprojzrram  DCROUT 


This  subprogram  Is  called  by  NSWIT.  Inputs  to  DCROUT  are 
the  six-by-six  sensitivity  matrix  and  the  final  co-state  error 
vector.  The  sensitivity  matrix  is  Inverted  and  the  error  vector 
is  premultiplied  by  this  matrix.  The  output  of  this  subprogram 
is  the  vector  quantity  by  which  the  initial  co-state  should  be 
changed. 

Subprogram  PRNT 

This  subprogram  is  called  by  NEWIT,  Values  of  the  inde- 
pendent and  dependent  variables  of  the  iterator  are  printed  on 
each  iteration. 

Subprogram  TRAJB 

This  subprogram  is  called  by  NEWIT.  A trajectory  is  calcu- 
lated with  the  initial  conditions  as  determined  by  NEWIT,  This 

i 

subprogram  also  calculates  the  error  in  the  final  co-state 
conditions. 

Subprogram  HK4 

This  subprogram  is  called  by  NEWIT  and  the  MAIN  PROGRAM. 
This  subprogram  is  basically  the  4^^  order  Runge-Kutta  integra- 
tor of  the  IBM  Scientific  Subrountlne  Package  without  the 
accuracy  checks.  The  SSP  integrator  was  changed  slightly  to 
allow  the  terminating  energy  condition  to  be  satisfied. 
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Subprogram  FIGNC 

This  subprogram  Is  called  by  RK4.  This  subprogram  eval- 
uates the  right  hand  side  of  the  system  of  first  order  differen- 
tial equations  of  the  state  and  co-state.  The  effects  of  the 
Earth *s  orbit  about  the  Sun  Is  considered  in  this  subprogram. 

The  control  variables  ( p and  ^ ) are  also  calculated. 

Subprogram  OUTNC 

This  subprogram  Is  called  by  HK4,  The  limiting  energy 
condition  Is  checked  In  this  subprogram.  If  the  terminating 
energy  condition  Is  exceeded,  a new  time  step  for  the  HK4  sub- 
program Is  determined  such  that  the  terminating  energy  condition 
Is  satisfied.  All  output  Is  changed  from  reference  orbit 
canonical  units  to  MKS  units. 

Subprogram  OPCOV 

This  subprogram  is  called  by  the  MAIN  PROGHAM.  A summary 
of  the  converged  trajectory  Is  printed, 

A flow  diagram  of  the  subprograms  of  SOLCO  Is  shown  below 
In  figure  D-1, 


5 


SCLCC 


Id* 

> 

G 

o:  j 

H 

• 

K < 

X 

U CO 

o 

Uj 

-x 

GU3 

o 

< < 

O IM 

> 

m 

a j 

» w 

ex 

H 

HC 

OH  U. 

m 

a 

z 

(/J 

iL  in  o 

ct 

H 

J 

• • I 

< 

o 

cr 

< < 

ox 

m 

a 

X 

H < • 

CO 

-) 

mo  a 

m 

< 

s 

KZ 

X • « 

H 

a 

o 

H m G 

H 

z 

c 

• s z 

X 

G 

< • • 

s 

G 

in 

UJ 

HU. 

< 

M 

iH 

im  H (\j 

u 

a 

»-u. 

X • 

z 

> 

< 

HP)  • tl 

< 

H 

LL 

ou. 

• u z 

a 

U 

z u 

z 

u 

u 

H 

UJ 

♦ CO  CL  X • 

<L 

O 

U CJ 

*-a  ♦ Gx 

cu 

lU 

UJ 

z < 

w <a 

cr 

X 

“5 

i-i  a 

< ►-  ♦ 

H 

< 

CL  (0 

IX  H o ^ • 

G 

a 

m m u.  ^ o 

a 

H 

q:< 

Lyco  •0*^«K 

M 

o 

a 

UJ  • O w-  • ^ 

H 

CL 

H 

G 

U.U- 

H »-•  O > a CL 

Z 

< 

u 

z u.  H a a a 

z. 

a 

LL 

a 

X 

*-•  Q.  m • UJ  z • 

I-M 

> 

a — 

<IUj 

• a Q CO  Uj 

X 

• p) 

H 

QC  in 

• Q *0 

D 

a 

m o 

a 

o<. 

' < CL  ^ • < 

j:  h 

t\j 

“^C 

o 

cx 

N a m i/j  (\j  CL  a 

< 3 

< 

ex  a 

1 • » X m H 

u. 

LU 

cr  • 

G 

a 

O m m • w H-  • 

(/)Z 

r' 

X 

HvO 

UJ 

CJ 

• O X > O')  o 

1-  •- 

H 

• ^ 

o 

oa 

X • o m * z z 

z 

a 

« cu 

LC 

z< 

1 OJ  Z • H 

< s 

h- 

in 

ZH 

CL 

<u<  'Si  IPX  G 

1-  < 

m 

* M 

> 

-Jin 

w • • • w < Q 

W)  £X 

U 

> 

H cr 

Z 

JUj 

TO-,.^^HX  • 

20 

o 

G 

z s 

u 

a 

* O X O X ^ O 

cu 

O 

G-* 

o 

ct  > 

J .o^i^o:  zz 

u a: 

in 

uo 

Her 

<ozxa\o 

a 

• <t 

^ U' 

m 

z< 

m < < N s z 

U-Q 

GCJO 

s 

WC^ 

X 

OH 

a WK  itcu. 

oc 

omo 

HO' 

H 

UlU 

<CD-<II.U 

(/I 

z 

2 

»-4  «! 

z 

t-  S N N N V _l 

in 

ax* 

so 

U. 

a;< 

•-*  4 

O UJ 

Z < >X) 

S 

^ m UJ 

o, 

i-j 

(J2Z2ZZZ 

ZI 

< 

II  z • 

na 

-oouaoo: 

►-I- 

H-  UJ 

a. 

a — 

> 

-IS  5 ? S S LI 

t- 

J z H 

IS 

•-  _JZ 

X 

in  UJ 

asiksssK 

h-  > 

G 

c 

_Jw, 

< 

^ X 

s o u c U C X 

U.C 

< X a 

a 

W4a 

X 

H 

“-OUUOOiiJ 

cn 

o^^  s 

a 

— u<-« 

X 

in 

u 

G. 

_1Z 

LJ 

CO 

xc 

-ID 

1 

I 

HU. 

<Q 

1- 

< 

• 

> 

d 

u 

H 

u 

U' 

— n 

> < 

X 

a.  cr 

G 

O H 

<. 

H 

• 

O G 

o 

Lj  *-*  < 

z 

X 

H 

< G G ^ 

G 

Ct  G G H 

u 

H G iL  a 

* 

X CO  c 

o 

G l/)  CO 

z 

^ CO  m m 

o 

a m o X 

w 

-OO-^H 

X 

HOG* 

* 

a G in  ioct 

X 

C CO  z c 

G 

H •-'X 

X 

X Z U O 

O G z m CO 

• 

a xm 

X 

> z G 

►- 

tl  UJ  U Z G 

a < ^ o < 

z 

H X H *-  > 

* 

U)  <►- 

> 

*-  X N < G 

11  u 

X ^ (1  m 

m H 

tr  X mo 

o 

m o H a 

a.  u H « oj 
«q:  a > 
a a uj  z 


H O 

-* 

XC 

OX  ^cu 

ro 

C H G 

ox  *0 

o 

XXX 

*-a  H o 

o 

oo  a 

• W O -4 

m p^  CM  • • 

UVJUU 


uuu  wuu 


sC  -St  • • • 

O w ^ u.  VO  (Vi 
II  a-)w  w 

a u.  X uj  G uj 
I-*  K J ^ H H ( 

^1-,  Jw«  ^1 

u)  a < u.  cr  c cri 
•-  5 U so  a < 


<<■<<< 

C.  2 X 2 S 5 
ccr  u:  X cr  a c 

KUUOUU-: 
cniL  u.  iL  u.  u.  ui 


^ CvJ  ^ CJ  if. 

ooooo 
o o o oo 


INPNCL 


WRITE( 6.1 C07) 
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